Ab Initio Calculation of Spin Gap Behavior in CaV 4 Og 
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Second neighbor dominated exchange coupling in CaV409 has been obtained from ab initio density 
functional (DF) calculations. A DF-based self-consistent atomic deformation model reveals that 
the nearest neighbor coupling is small due to strong cancellation among the various superexchange 
processes. Exact diagonalization of the predicted Heisenberg model yields spin-gap behavior in good 
agreement with experiment. The model is refined by fitting to the experimental susceptibility. The 
resulting model agrees very well with the experimental susceptibility and triplet dispersion. 
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CaV40g was the first two-dimensional system observed 
to enter a low-temperature quantum-disordered phase 
with a spin gap A w 110K. The gap was first apparent 
in its susceptibility, which vanishes at low temperatures 
as x(T — » 0) ~ exp(— A/fcT) 0, and was observed di- 
rectly in the dispersion of triplet spin excitations (^q) 
measured by neutron scattering 0] . This unexpected be- 
havior has stimulated considerable theoretical study of 
the exchange couplings between S=^ spins on the V lat- 
tice using Heisenberg models [||-|| . 

CaVjOg is a layered compound — the interlayer dis- 
tance is sufficiently large to make interlayer V-V coupling 
negligible. Within a layer, the V atoms form a i-depleted 
square lattice shown as the circles in Fig. [j] The 
lattice was originally viewed as an array of square "pla- 
quettes" of V ions (e.g., 1-2-3-4 in Fig. [l]) tending to- 
ward singlet formation since isolated plaquettes have a 
singlet ground state. Examination of the structure how- 
ever suggests intra- and inter-plaquette nearest neighbor 
V-V coupling should be similar, so the limit of isolated 
plaquettes is not realistic. 

Self-consistent electronic structure work || identified 
the V + spin orbital as d xy , which implied that it was 
a larger square of V ions, the "metaplaquette," where 
singlet formation arises. Fitting Heisenberg Hamiltoni- 
ans to the measured dispersion of the triplet excitations 
confirmed that the dominant second neighbor exchange 
coupling is crucial to account for the shape of Qq || . 

The complete Heisenberg Hamiltonian for CaV^g has 
four different coupling constants: nearest-neighbor (nn) 
and next-nearest-neighbor (nnn) couplings and, for each 
of these, intra- and inter-plaquette couplings. In notation 
of Gelfand et al. || , the Hamiltonian is given by 

H = Ji Si ■ Sj + J[ Sj • Sj 

nn nn' 

+ J 2 ^S i -S J + J^Si-S,-, (1) 

nnn nnn' 

where Si denotes the spin ^ operator in site i. The nn 




FIG. 1. Couplings in CaV4 0t). The circles represent V 
atoms, and the lines between V's show the couplings. The 
numbers label the sites used in the LSDA calculations to de- 
termine the couplings. Line thicknesses are proportional to 
the best set of couplings found by fitting to the experimen- 
tal susceptibility. The strongest coupling, J' 2 , is shown as the 
thick dot-dashed lines, forming metaplaquettes, e.g., 1-6-3-8. 

sums run over nearest-neighbor bonds and the nnn sums 
run over next-nearest-neighbor bonds. Unprimed sums 
connect V's in the same plaquette, while primed sums 
connect V's in different plaquettes. The four couplings 
are drawn in different line styles in Fig. |l|. 

In this Letter we show that the spin gap behavior 
of CaV40g, even considering its complex structure with 
eight very low symmetry V 4+ ions in the primitive cell, 
can be calculated in ab initio fashion. Our work has 
three separate aspects. 1) Local spin density approxi- 
mation (LSDA) calculations are used to obtain energies 
for various magnetic configurations. The resultant ex- 
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TABLE I. Magnetic configurations of the eight V ions in 
the primitive cell for the states used to determine the ex- 
change constants from LSDA. Most configurations are defined 
in the text. V ions are numbered as in Fig. 1. The final col- 
umn shows the relative LSDA energies. 
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change interactions are obtained by fitting these ener- 
gies to the mean-field Heisenberg model as described be- 
low. 2) An approximate but physically motivated local 
orbital method called the self-consistent atomic defor- 
mation (SCAD) method j|] is used to provide explicit 
local orbitals, eigenvalues, and hopping integrals for cal- 
culating the exchange interactions from perturbation the- 
ory This method reveals that the nn interactions are 
not intrinsically small, but the net value of the superex- 
change coupling is small due to cancellations among var- 
ious fourth-order processes. It also indicates that direct 
V-V exchange coupling is important. 3) The Heisenberg 
Hamiltonian is solved using exact diagonalization tech- 
niques on finite periodic clusters. Spin gap behavior is 
obtained, and x(T) is similar to the data. The Heisen- 
berg couplings are refined by fitting to x(T)- The re- 
sulting Hamiltonian agrees well with x(^) arl d with the 
triplet dispersion determined from neutron scattering. 

The LSDA calculations of the energy for various mag- 
netic configurations were more precise extensions of pre- 
vious work on CaV^g p,^0[. The magnetism of the V 
ion is found to be robust, allowing us to break the spin 
symmetry in any manner we choose and obtain the en- 
ergy from a self-consistent calculation. The symmetry 
of the non-magnetic state is initially broken as desired 
by applying the necessary local magnetic fields to the V 
ions. The seven configurations we have chosen include 
the ferromagnetic (FM) state, one ferrimagnetic (FiM) 
state, and five antiferromagnetic (AF) states with zero 
net spin. These AF states include the Neel state, a state 
in which FM plaquettes are antialigned (FMPL), and a 
state in which the metaplaquettes are aligned antiferro- 
magnetically (AFMP). The configurations, given explic- 
itly in Table |, were chosen either because of their phys- 
ical relevance (AFMP was anticipated to be lowest in 
energy, as found) or computational considerations such 
as retaining inversion symmetry. 

The resulting energies were fit to the mean-field 
Heisenberg model, which contains simply the Sf or Ising 
terms of the full Hamiltonian (|l|), to determine the four 
coupling constants. The six energy differences lead to six 



TABLE II. Values for the four couplings (in meV). The 
LSDA values are derived from the energies in Table I. The 
SCAD results are derived from the local orbital method, and 
the Fit results come from fitting the experimental suscepti- 
bility. Both are described later in the paper. Also shown are 
the couplings deduced from neutron scattering data (Jm. 
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conditions on the four Js, and a least-squares fit gives 
the values listed as LSDA in Table ||, each with a fitting 
uncertainty of about 1 meV. Since both nearest and next 
nearest couplings are AF in sign, there is a great deal of 
frustration in the magnetic system. The large value of 
J' 2 indicates that singlet formation on the metaplaquette 
is the driving force for the spin gap. 

To understand how these values of the exchange pa- 
rameters arise, we evaluate the fourth-order expressions 
for the exchange constants, using an approximate but 
parameter-free method based on the SCAD method. For 
each coupling constant in CaV40g, we focus on the rel- 
evant clusters for each coupling. The nnn interactions 
require a V2O cluster with two V ions (each with one rel- 
evant orbital) and one O in between. The nn exchange 
interactions require a V2O2 cluster. All three 2p orbitals 
in each O are relevant, since the low symmetry makes 
them non-degenerate and oriented in directions deter- 
mined not by symmetry but by electronic interactions. 

We neglect the Hubbard U and Hund's rule coupling on 
the O ions. In what follows, U is the V on-site repulsion, 
ey and e a are site energies of the V and a-th O orbitals, 
and ti a is the hopping amplitude between the i-th V and 
the a-th O orbital. Defining the energy denominators 
A Q = U + ey — e Q simplifies the expressions. 

The initial state has each O orbital doubly filled and 
each V with one electron. The perturbation theory is 
given by three fourth-order terms and the direct second- 
order V-V term: 



J = h + h + h + jd 




In the nnn case, a and (3 sum over the three orbitals in 
the single oxygen atom. In the nn case, a and (3 sum over 
the six orbitals in both oxygen atoms. The first three 
terms in (|^) can be categorized by their configurations 
after the second hop of the four-hop process: 1) One 
vanadium empty; 2) One oxygen orbital empty; 3) Two 
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oxygen orbitals half filled. The last term has an extra 
factor of two because it arises twice: the total spin singlet 
case is reduced in energy and the total spin triplet is 
increased by the same amount. The latter picks up a 
minus sign due to electron exchange. 

This expression is evaluated with the SCAD model, 
which expresses the total density n(r) as a sum over lo- 
calized densities |(^(r — R;)| 2 centered at the atomic 
sites Ri |)|. The orbitals <fia^ are solutions to atom- 
centered one-electron Hamiltonians Hi for each site. The 
potentials in Hi are determined self-consistently from the 
expression for the functional derivative of the total en- 
ergy. It includes a local approximation for exchange and 
correlation energy Jll| l and the Thomas-Fermi function 
for kinetic energy of overlapping densities. 

Each V ion has the lowest of its five 3d levels occupied 
by a single electron, giving the V 4+ , 2 ~ ionic descrip- 
tion. U ~ 3.5 eV was computed by minimizing the SCAD 
energy subject to the constraint that one V ion has its 
charge increased by unity. The electron comes mainly 
from the other V ions with only a minor portion coming 
from the nearby O ions. 

The matrix elements, t%j — {ipi\H\ipj) require the full 
Hamiltonian H and orthogonalized orbitals ip. The ^'s 
are obtained from the SCAD orbitals using Lowdin's 
method JT2| , and H is determined from the site centered 
SCAD Hamiltonians by removing the kinetic energy over- 
lap contributions from the latter. This gives expressions 
for H that differ in the site selected for spherical har- 
monic expansion of the potential. We find the two possi- 
bilities, and tji, may differ by ^20%, which leads to a 
much larger uncertainty in the fourth-order J's. Since the 
vanadium sites of a given pair of V ions are equivalent by 
symmetry, the direct interaction, jd, has no such uncer- 
tainty. To be consistent with the direct interaction cal- 
culation, we use the vanadium-site-expanded potentials 
for evaluating matrix elements between oxygen- vanadium 
pairs. The net values obtained (labelled SCAD in Table 
|J) agree rather well with those derived from LSDA en- 
ergies for Ji, J 2 , and J' 2 - The close agreement may be 
fortuitous in view of the uncertainties mentioned above 
and the approximations inherent in the SCAD method. 
Nevertheless, we believe certain qualitative features of 
the SCAD results are real: 1) The values for J\ and J[ 
result primarily from jd, with relatively small contribu- 
tions from fourth-order terms due to cancellation within 
ji and between j 2 and j'3. 2) The value for J 2 , the largest 
coupling, is dominated by a single term in ji, resulting 
from V overlap with the middle O 2p level. 

For each set of four coupling constants, we calculated 
the uniform susceptibility of the Hamiltonian (|l|) on pe- 
riodic 20-spin clusters. The susceptibility is given by: 
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FIG. 2. Uniform magnetic susceptibilities calculated by 
exact diagonalization of a 20-spin cluster. The theoretical 
curves using the coupling constants from Table [n| are shown 
as lines, while the circles show the experimental susceptibility 
of Taniguchi, et at The theoretical fit to the susceptibility 
is the solid curve that lies over the experimental points. 



where n is the number of V atoms per gram and N is 
the number of sites in the cluster. We take g — 1.67 for 
all plots. This was determined from the fit to the exper- 
imental magnetic susceptibility data described below. 

To evaluate (|^), we calculate all eigenvalues of the 
Hamiltonian — eigenvectors are not required. We block- 
diagonalize the Hamiltonian with all possible symme- 
tries: translations, rotations, S, and S z [l3|| . The blocks 
are left with no degeneracies, so the eigenvalues are 
calculated very efficiently using the Lanczos algorithm 
with no rcorthogonalization developed by Cullum and 
Willoughby ]l4| . This allows x to be calculated exactly 
at all temperatures using one Lanczos run for each sym- 
metry sector. The Hamiltonian for the 20-spin cluster has 
blocks as large as 36950. Within each block at least the 
400 lowest and highest eigenvalues are calculated, and 
an analytic density of states is assumed for the middle 
eigenvalues. This technique will be described elsewhere. 

The susceptibility of the full Hamiltonian (|l|) calcu- 
lated with each set of coupling constants in Table || 
is shown in Fig. 2. The experimental susceptibility of 
Taniguchi, et al. [1] is shown for comparison. All curves 
exhibit a spin gap, as evidenced by their low temperature 
behavior, \(T —> 0) ~ exp(— A/kT), where A is the gap. 
Both the LSDA and SCAD approaches overestimate the 
gap, indicating that the calculated coupling constants are 
too large. The coupling constants deduced from neutron 
scattering are also shown [f|||]. 

Also shown in Fig. || is a curve generated using the 
coupling constants obtained from a least-squares fit of 
the susceptibility to the experimental results. In the fit- 
ting procedure, we allow the g-value in eq. (|§|) and all 
four J's to vary. At the best fit, we obtain the cou- 
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FIG. 3. The triplet dispersion £Iq in CaV^Og, calculated 
from the fifth-order metaplaquette series expansion of Wei- 
hong, Oitmaa, and Hamer ||. The circles are the neutron 
scattering data of Kodama, et al. || . The ab initio couplings 
have been rescaled so their minimum gaps match the experi- 
mental minimum. The solid line shows the dispersion of the 
(unrescaled) couplings determined by fitting the experimental 
susceptibility. 

pling constants listed as "Fit" in Table [n] and shown as 
the line thicknesses in Fig. 1. We find g = 1.67, which 
is smaller than the g-value indicated by ESR measure- 
ments pi . Near the minimum, the fitting function is 
quadratic. The eigenvalues of the Hessian (scaled by 
an arbitrary constant) are 1, 0.046, 0.013, and 0.00039. 
The smallness of the last eigenvalue indicates that in 
the 6{Ji, J{, J 2 , 4} = {0.09,-0.57,0.81,-0.09} direc- 
tion from the minimum, the least-squares fit is very soft. 

The 20-spin cluster is sufficiently large compared with 
the correlation length to describe the infinite system ac- 
curately. The minimum triplet gap hardly varies be- 
tween 20 and 32-spin clusters: A20 = 9.92 meV while 
A 32 = 10.02 meV for the Fit Hamiltonian. 

Fig. |^ shows the triplet dispersion J7q of the LSDA, 
SCAD, and susceptibility-fit coupling constants calcu- 
lated with the expansion in Ref. ||. Since the LSDA 
and SCAD coupling constants overestimate the gap, we 
rescaled their J's by 0.58 and 0.65, respectively. Both the 
Fit and rescaled LSDA J1q agree with the neutron scat- 
tering data reasonably well; in particular, they correctly 
have minima at Q = (0,0). 

To conclude, we have shown that the quantum- 
disordered phase in CaV^g can be predicted in ab ini- 
tio fashion. We calculated the coupling constants of 
the Heisenberg Hamiltonian for CaX^Og in two very 
different first-principles approaches. In both methods, 
the strongest coupling is found between next-nearest- 
neighbor V atoms on metaplaquettes — the weak coupling 
between nearest-neighbor V's results from the cancella- 
tion among superexchange processes. The uniform mag- 
netic susceptibility for each set of coupling constants is 
calculated using a novel finite-temperature exact diago- 



nalization technique, which shows the Hamiltonians de- 
termined from both ab initio approaches have quantum- 
disordered phases. The Hamiltonian that best fits the ex- 
perimental susceptibility is calculated, and the agreement 
is remarkable. Finally the triplet dispersion of the ab ini- 
tio and best susceptibility-fit Hamiltonians are shown to 
agree well with the neutron scattering data. 
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